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Abstract 

I describe the Time-Dependent Superfluid Local Density Approximation, which is an adiabatic 
extension of the Density Functional Theory to superfluid Fermi systems and their real-time 
dynamics. This new theoretical framework has been applied to describe a number of phenomena 
in cold atomic gases and nuclear collective motion: excitation of the Higgs modes in strongly 
interacting Fermi superfluids, generation of quantized vortices, crossing and reconnection of 
vortices, excitation of the superfiow at velocities above the critical velocity, excitation of quantum 
shock waves and domain walls in the collisions of superfluid atomic clouds, excitation of collective 
states in nuclei. 
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1 Density Functional Theory 

In a remarkable theorem proven almost fifty years ago, Kohn and collaborators 
'Hill' established that there is an one-to-one map between the ground state wave 



function of an interacting multi-electron system, the number density, and the ex- 
ternal Coulomb one-body potential created by the (static) nuclei. Even though 
the theorem was formulated in terms of electrons and Coulomb interaction among 
them and with nuclei, the proof never relies on these specific aspects and it ap- 
plies to any non-relativistic fermion system and as such it is widely referred to 
as the Density Functional Theory (DFT). Subsequently this theorem has been 
extended to apply to more general situations (|4|[5|) and since then it has became 
a workhorse in chemistry and condensed matter studies, due to the tremendous 
mathematical simplification achieved by replacing the many-body Schrodinger 
equation with a system of non-linear and coupled 3D partial differential equa- 
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tions, formally equivalent to a meanfield treatment of the electron systems. The 

mam difficulty in applying DFT resides in the fact that the generation of the cor- 
responding Energy Density Functional (EDF) is more of an art than a science, 
as the Hohenberg-Kohn (V) and Kohn-Sham Q theorems do not provide for 
an algorithm to derive the EDF from the many-body Schrodinger equation. One 
very important generalization of the original DFT approach was the development 
of the time-dependent version of DFT ([6|[7|[8]), which allows in principle the re- 
placement of the time-dependent many-body Schrodinger equation with a set of 
time-dependent coupled 3D partial differential equations, formally equivalent to a 
time-dependent meanfield approach. Thus "Time-dependent density functional 
theory (TDDFT) can be viewed as an exact reformulation of time-dependent 
quantum mechanics, where the fundamental variable is no longer the many-body 
wave- function but the density." Q. Since its formulation TDDFT has been ap- 
plied mainly in atomic and molecular calculations for the study of the excited 
states of many-electron systems and by now it has achieved a very sophisticated 
level of development as amply highlighted in recent monographs ([7||8|). While 
in principle DFT as initially formulated can be used to describe any interacting 
fermion system, the application to superfluid fermionic systems is clearly going 
to encounter some challenges. If one were to use only the density of the many 
fermion system it would be impossible to distinguish between a superfluid and 
a normal system. One can easily envision a situation when a superfluid system 
is stirred energetically enough that in some regions or even globally the system 
can undergo a transition to a normal state. The need to develop a version of the 
DFT suitable for the description of superfluid systems was recognized quite some 



time ago by Oliveira, Gross and Kohn (10). The version of DFT suggested by 
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these authors however lacked the great simphcity of the Local Density Approxi- 
mation (LDA) of Kohn and Sham (2), namely of being formally equivalent to a 
local implementation of the meanfield approximation. Similarly to the approach 



of Bardeen, Cooper and Schrieffer (11) Gross et al (10 ,12 ,13) introduced the 
anomalous density to describe the order parameter in the superfluid phase. The 
anomalous density however has an ultraviolet divergence if the pairing field is 



a local potential (14 15 ,16 17) and a meaningful theoretical framework can be 
developed only if one introduces and appropriate regularization and renormaliza- 
tion procedures, similar to those routinely used in quantum field theories. Gross 



et al (10, 12, 13) eschewed these issues by considering a non-local version of the 
DFT in the case of superfluid systems, with a non-local pairing field. 

In the case of Bose superfluids, "two phenomenological theories explain almost 



all experiments to date" (18), the two-fluid hydrodynamics and the "comple- 
mentary view, provided by Fritz London, Lars Onsager, and Richard Feynman, 



(that) treats the superlfuid as a macroscopic quantum state" (18). Landau (19) 
identified the excitations which he called rotons with motion characterized by 
non-vanishing superfluid velocity circulation (V x -u 7^ 0) and he did not envision 
the existence of quantized vortices. When a superfluid is brought into rotation, 



quantum vortices, predicted by Onsager and Feynman (20,21), are formed and 
the two-fluid hydrodynamics is unable to describe their dynamic generation. Even 



though Landau described his approach initially as a quantum theory (19), only in 
its classical incarnation has ever been used in practice, and the two-fluid hydro- 



dynamics (23,24,19,22) is "essentially thermodynamics" (18). It is ultimately a 
phenomenological classical approach in which Planck's constant never explicitly 



enters. Later on, vortex quantization was imposed by hand when needed (22), in 
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a manner similar to the Bohr quantization rule of the hydrogen atom, following 



Onsager's and Feyman's quantization conditions (20,21). 



In the case of dilute Bose systems the Gross-Pitaevskii (25,26) equation can 
be derived 

ih'^{f,t) = - — ^{r,t)+g\^{f,t)\'^^{r,t) + Ve^t{r,t)^{r,t), (1) 
2m 

where g' > is a coupling constant for the boson-boson interaction and Vext i^, t) 
is an external potential. This equation can be used to describe the real-time 
evolution of a bosonic superfluid at very low temperatures, and in particular the 
vortex generation. A similar description did not exist for a fermionic superfluid 
until recently. Many phenomena in cold atom physics, nuclear physics and neu- 
tron star crust require a real-time dynamical approach, often beyond the linear 
response regime and for time scales when the role of collisions can be neglected. 
The description of such phenomena lead the need for a theoretical framework of 
the real-time dynamics of a fermion superfluid arises. 

A local extension of DFT to superfluid systems in the spirit of Kohn-Sham LDA 



has been formulated in Refs. (16,17,27) and was dubbed the Superfluid Local 
Density Approximation (SLDA) as being a natural extension of the Kohn-Sham 
LDA to superfluid systems. SLDA and its time-dependent extension TDSLDA 
was applied to a range of static and time-dependent situations in nuclear and 
cold atom systems ((27|[29|[30}|3ll[32|[33}[3§[35^ In its present formulation 

SLDA has been implemented only for phenomena at zero temperature, and its 
extension to time-dependent processes is valid only as long as collisions, leading 
to entropy production and thermalization of the system can be ignored. In this 
respect TDSLDA is similar in spirit to the Landau's Fermi liquid theory, which 
describes the propagation of the high-frequency zero-sound at not very long time 
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scales, as opposed to a kinetic approach necessary to describe the propagation 

of the relatively lower frequency first-sound at much longer time scales. In the 
case of the zero-sound the local Fermi distribution does not equilibrate, while 
the opposite is true in the case of the first-sound. The high-frequency and low- 
frequency are determined when compared to the local relaxation rate, governed 
by the collision integral in a Boltzmann description of such system. In the small 
amplitude limit the TDSLDA equations in the frequency representation become 
formally equivalent to the Landau's Fermi liquid theory or the linear response 
theory in the presence of pairing correlations (or their absence, if the pairing field 
vanishes) . 

2 Unitary Fermi Gas 

In 1999 Bertsch introduced a hypothetical system which later became known as 



the Unitary Fermi Gas (UFG) (38). At the time this was a pure theoretical 



model of very dilute neutron matter, which subsequently became an object of 



intense study both theoretically and experimentally in the cold atom physics (39 
|40||4ll|42||43|[28j). The UFG has properties surprisingly close to those of realistic 
dilute neutron matter (44:). A UFG is a system of spin-1/2 fermions interacting 
only in the s-wave with an infinite scattering length and a zero effective range. 
The only Trivial dimensional arguments show that the energy of a homogeneous 
UFG is a function of only the Planck's constant h, the fermion mass m, the 
volume of the system V and the number of fermions N. The only quantity with 
dimension of energy one can form out of these constants is x SepN/b, where 
is a dimensionless constant now called the Bertsch parameter, ep = h^kp/lm is 
the Fermi energy of a free Fermi gas of the same density, and N/V = kp/Sn'^. 
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Since such an interaction is attractive ^ < 1, but in 1999 it was not clear whether 

also > 0. For = 1 this is the energy of a free Fermi gas, and if ^ < 

the system would collapse. Accurate quantum Monte Carlo (QMC) calculations 

^ = 0.372(5) ([45}|46||47||48]) and experimental measurements ^ = 0.376(4) (|49) 



have converged to basically the same value. Consequently a UFG is a gas, but 
somewhat surprisingly is also a superfluid, with one of the largest relative known 
pairing gaps in any fermion system A ~ O.bep. Superfludity in a UFG has been 



confirmed experimentally (50) by directly putting in evidence the formation of 



the Abrikosov lattice of quantum vortices (51), when such a system is brought 
into rotation by stirring it with laser beams in an atomic trap. 

The (TD)SLDA equations of motion are derived using an appropriately de- 
fined action integral, which is a straightforward generalization of the Kohn-Sham 
approach ^ to superfluid systems and time-dependent phenomena. This kind 
of extension of the Kohn-Sham approach to time-dependent phenomena is ap- 
propriately referred to as the Adiabatic Local Density Approximation (ALDA). 
If the number of spin-up and spin-down particles are equal, there is no spin-orbit 
interaction or velocity coupling, and if spin degrees of freedom are not excited, 
this action integral has the form: 

S = in J dtcfr {< (f, t)dtUn{f, t) + < (f, t)dtVn{r, t)} 
+ / dtd^rl £ 



n(f,t),T(r,t),J(f,t),z.(f,t)J +U{r,t)n{r,t)^, (2) 

where Un{r,t) and Vn{r,t) are quasi-particle wave functions (qspwf(s)), and S 
is the energy density, which depends on the number n{f,t) = 2J2n\'^nir,t)\'^, 
kinetic r(f,t) = 2 ^„ | Vu„(f, t)p, current j{f,t) = 2 Im[fn 2„ <(f, t)Vt;„(f, t)], 
and anomalous density ^{r^t) = X^n ^n(^) 0^n('^) i) respectively. U{f,t) is some 
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external potential in which the system might reside. Sometimes it is convenient 

to consider as well and external pairing field acting on the system, by adding to 
the above action integral a term Aext{r,t)h'*{f,t) + A*^^{f,t)h'{f,t). 

The energy density £, as usual in DFT, is universal, in the sense that its form 
is independent of the external potential U{f, t). This is the reason why £ should 
satisfy general symmetry principles: translational and rotational invariance, par- 
ity, local gauge and local Galilean invar iance/covarince, and the theory should be 
renormalizable as well. As was mentioned above, in the case of superfluid systems 
the anomalous density and the kinetic energy density are both ultraviolet diver- 



gent quantities (14,15,16,17) and in order to remove these divergencies from the 
formalism both the kinetic and anomalous densities should enter the formalism 
in a unique combination. The UFG is a quire remarkable physical system, as one 
can use also simple dimensional arguments in conjunction with the symmetry 
requirements and renormalizability of the theory to show that the energy density 
has a very simple and essentially unique expression (if simplicity of the theory is 
invoked) : 



n(r, t), T(f, t),j{f, t), v{r, t) 



a 



— T{r,t)-A{r,t)u*{r,t) 
2m 



10m mn[r,t) mn{r,t) 

The first term is the unique combination of the kinetic and anomalous densities 
required by the renormalizability of the theory; the second term is the only func- 
tion of the density alone allowed by dimensional arguments in the case of a UFG 
(since no other dimensionfull parameters are needed to describe this system) ; the 
third term is lowest gradient correction and its amplitude is likely small ('52'); 
and the fourth term is required by Galilean invariance ([29j). The last term van- 
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ishes in the ground state where there are no currents, and the gradient correction 

term is absent in homogeneous matter. The constants a, (3, A and 7 (which is not 



shown here, see Refs. (27 ,29) for details, but enters in the definition of the pair- 



ing field A(r, t) (29)) are all dimensionless and they fully determine the energy 
density functional for a UFG. The constants a, /3, 7 are determined by requiring 
that this energy density reproduce exactly the energy per particle, paring gap 
and effective mass obtained in accurate QMC calculations for a uniform system 



(17,45,46). The constant A is determined from reproducing the QMC results of 



an inhomogeneous UFG in an external trap ( 52 ) . The magnitude of the gradient 
corrections is relatively small A « —0.1 (and somewhat surprisingly negative) 
(52) and the effective mass is close to the bare mass value, as q ~ 1.1 (|27|[29}[52]). 
Consequently, the last two terms in Equation ([3]) represent somewhat small cor- 
rections when compared to the first two terms. Using this EDF (without gradient 
corrections, which were introduced later to further improve the agreement) one 



can predict now with impressive accuracy (27,29) the results of the QMC cal- 



culations of inhomogeneous systems (53,54) without any additional fitting and 
the agreement is always within the QMC errors. Even though more complicated 
EDF could be contemplated, which however will lack the simplicity of Equation 
([3]), there does not seem to be any need for them yet at the level of accuracy of 
the present QMC calculations. The absence of any dimensionfull scales, apart 
from the average inter-particle distance n^^^^ « T^/kp, allows the determination 
of the structure of the EDF for a UFG essentially uniquely, and sets this strongly 
interacting system apart from other many-body systems. 

The equations for the qpwfs Un{r, t), Vn{r, t) have the time-dependent Bogoliubov- 
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de Gennes (55) form by design: 



ihdtUnif,t) 
^ ihdtVn{f,t) J 



^ h{f,t) + Uir,t) A{f,t) \ I Unir,t)^ 



■ (4) 



\^ A*{r,t) -h*{f,t)-U{f,t) J \vn{r,t) J 
The operators here are defined as corresponding variational derivatives as usuah 
h{f,t) = 58 / 8n{r,t) and A(r, t) = 6£/5v*{r,t). The fact that the pairing po- 
tential A(r, t) is a local operator, and not an integral one, is the reason why a 
Kohn-Sham type of implementation of DFT for superfluid systems displays ultra- 



violet divergencies in the kinetic and anomalous densities (14,15,16,17), which 
require the implementation of a regularization and renormalization procedures 
described in Refs. ( [T6|[T7|[27|[29| ) . Numerically, the solutions of these equations is 
a formidable problem, since one has to solve an infinite system of time-dependent 
3D nonlinear coupled partial differential equations, and this is feasible only on 



modern supercomputers ( 56 ) . 



3 Nuclear Systems 

During the last decade it was realized that the traditional nuclear meanfield-like 
approaches used for decades are nothing else but a disguised form of DFT and 
the corresponding terminology has entered the theoretical nuclear physics field as 
well. The construction of an accurate nuclear EDF is an ongoing endeavor and 
several approaches are used. A number of nuclear theorists hope to handle the 
strong nuclear interactions within a many-body perturbation theory, using the 
modern chiral perturbation theory nucleon interactions to construct the ground 
state energy of nuclei ( [57|58p9l ), from which one can derive accurate enough local 
nuclear EDFs. Typically, however, various authors prefer a more phenomenolog- 
ical approach, using local EDFs which depend on proton and neutron number 
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densities, spin and kinetic energy densities (and some of their derivatives), cor- 
responding current densities, and proton and neutron anomalous densities as 
well (|60|[6l||62|[63l|64]). These nuclear EDFs should satisfy the usual symme- 
tries: translational and rotational symmetry; isospin symmetry which is broken 
explicitly by proton-neutron mass difference. Coulomb interaction and charge 
symmetry breaking forces (this one being routinely neglected in practice) ; parity, 
gauge, and Galilean invariance/covariance, and renormalizabilty (often treated in 
a rather naive manner). Unfortunately, a number of these general requirements 
are often sacrificed in practice in the name of reaching agreement with experi- 



ment, see Ref. (65) for a discussion of some these issues and of potential ways 
to further develop TDDFT. In nuclear systems there is a strong spin-orbit cou- 
pling, and both neutron and proton systems can become superfluid. In principle 
it is also possible that proton-neutron Cooper pairs are formed in some instances, 
though the experimental evidence is ambiguous at this time. It is possible as well 
that neutrons form pairs in p— and /-waves (coupled by spin-orbit interaction) as 
well at higher densities in the neutron stars and the corresponding pairing field 
has a rather complicated spin-orbital structure. However, the basic difference 
of nuclei from UFG is the the presence of two types of particles (protons and 
neutrons), spin-orbit coupling and consequently the dependence of the EDF on 
several kinds of densities. However, formally the corresponding equations for the 
qpwfs have a similar structure. 
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4 Excitation of various collective modes 



4.1 The Higgs mode 

The Higgs mode in fermionic superfluids is perhaps one of the most intriguing 
cohective excitation modes of such a system, as its characteristics defy our usual 
concepts about collective motion. If one were trying to stretch or compress adia- 
batically any system, the work performed would be interpreted as the potential 
energy corresponding to that collective degree of freedom. One would have to add 
to that a corresponding collective inertia and with the emerging collective Hamil- 
tonian one would be able to predict rather accurately the oscillations (sometimes 
of both small and large amplitude) around the ground state equilibrium configu- 
ration. In the case of a fermionic superfluid one characterizes the work performed 
by an adiabatic "external" pairing field Aexti'f^,t) with a potential energy corre- 
sponding to this collective mode having the shape of a "Mexican hat," e.g. the 
energy of a superconductor just below the critical temperature in the Landau- 
Ginsburg phenomenological model. The "Mexican hat" potential is ubiquitous 
in describing various symmetry breaking mechanisms. The Anderson-Bogoliubov 



sound waves, with a dispersion uj = ck (66,55), represent the Goldstone modes 



which appear as a result of the broken gauge symmetry in the case of a superfluid. 
The Goldstone modes are described by waves of the phase of the pairing field, 
which propagate along the bottom of the Mexican hat potential and require a 
very small energy, see Figure Q. The oscillations of the phase of the pairing 
field are accompanied by small oscillations of the number density as well. In the 
case of a UFG the sound velocity is given by c = \/^fikp / \/3m, which differs 
from the free Fermi gas only by the factor < 1. The Higgs mode on the other 
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hand corresponds to radial oscillations of the pairing field along the radius in 

the Mexican hat potential. Naively one would expect that if the pairing field is 

excited in the radial direction with a small amplitude, one would observe radial 

oscillations of the pairing field with a finite frequency (unlike sound modes, whose 

frequency tends to zero in the long wave-length limit), predicted to be exactly 



hu = 2Aq (66), where Aq is the the ground state value of the pairing field. It 



was realized later on however that this is not correct, and Volkov and Kogan (68) 
have shown that the oscillations of the pairing field couple with excited quasi- 
particles with energies above the "new" gap 2Aoo < 2Ao, and that leads for large 
times to oscillation of the pairing field A(i) = A^o + Asm{2Aoot + Aoot. 
In this case the amplitude decreases very slowly cx l/\/t evolves towards a new 
equilibrium state with a smaller pairing gap Aqo < Aq and a certain fraction of 
excited quasi-particles. 

These results ( [68| ) were obtained for Fermi superfluids in the weak coupling 
limit, when the oscillations of the pairing field lead to no changes in the quasi- 
particle self-energy. In a UFG however any change in the pairing field induces a 
large change of the quasi-particle self-energy, since 

In the case of the Higgs mode however the number density stays constant n(r, t) = 
riQ, where no is the ground state value of the number density, and the oscillations 
have an infinite wavelength in a perfectly homogeneous system. Even though 
the TD-SLDA equations have a more complex structure, the dynamics of these 
modes (32 ) is very similar to that predicted for the weak coupling BCS superfiuids 
(68), see panel c of Figure [2] ([32j). However, if one excites a radial mode with 
larger amplitude, one observes a qualitatively different behavior, see panels a 
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and b of Figure [2j Surprisingly, if the system is taken almost adiabatically out 

of its equilibrium value Aq towards the tip of the Mexican hat potential and 

then released, the oscillations of the pairing field resembles a soliton train in 

time. Notice that most of the time the system illustrated in panel a of Figure 

[2] is essentially normal, and only for very brief periods of time the pairing gap 

attains values almost equal to the equilibrium one. The occupation probabilities 

of the single-particle states are very close to those of a normal system when 

A ^ Aq, but very different from their equilibrium values when the pairing gap 

attains values close to Aq, see upper two panels of Figure [3} Perhaps the most 

surprising feature is that the paring field while it oscillates, it never exceeds the 

ground state value Aq and it never swings on the other side of the minimum 

of the Mexican hat potential as one would have naively expected. The modes 

displayed in panels a and b of Figure [2] are truly nonlinear, non-equilibrium (see 

panels a and b Figure [s]) , and their frequency depends strongly on the oscillation 

amplitude, see the panels c and d of Figure |3] (32) . These are at the same time 

very slow modes, with frequencies well below the pairing gap H^Ih < 2Ao, but at 

the same time they are truly large amplitude collective modes, not only because 

of the size of the oscillation amplitude, but also because their excitation energy 

is equally large, see panel c in Figure |3] ([32]). In this respect these modes are 

somewhat similar to the very large amplitude oceanic waves, which also have 

very long wavelengths and even though carry a lot of energy, take a long time 

to dissipate into heat. The results of the somewhat related MIT experiment 



(69) suggests that the damping of these modes, due to the decay into other 
modes, is much smaller than one might naively expect at unitarity, related to 
the fact that the shear viscosity of the UFG is one of the lowest encountered 
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in Nature ([70|[7l}[72|[73]). A phenomenological Gross-Pitaevskii like description 

using Equation [l] for example, where ^'(r, t) would describe the order parameter 

(X A(r, t), would fail spectacularly. Since in this case there is no coordinate 

dependence (and in this case Aexti^^t) = ^extit)), the wave function ^{t) would 

simple have a phase evolution only and its magnitude would never change in 

time, unlike the behavior in Figure [2] ( |32[ ). The mechanism by which the Higgs 

modes dissipate their energy is still an open question, and so far they have no 

been put in evidence in experiments. A potential decay mechanism was suggested 

in Reference ( j74| ). In inhomogeneous systems the Higgs mode couples to density 



oscillations (32) and that is perhaps the best way to put them in evidence in cold 



atom systems. 

4.2 Nucleation and dynamics of quantized vortices 

Many more example of numerical experiments of vortex nucleation and dynamics 
discussed in this section can be found in Reference (33j) and especially on the 



webpage (75), where more than 4 hours of videos are available. The UFG is a very 
interesting system when it comes to its dynamics for one more reason, the Landau 
critical velocity, the flow velocity at which according to Landau a superfluid 
looses its superfluidity, is the largest known for any superfluid in appropriate 
units. There are two mechanisms relevant to our discussion which can lead the 
the loss of superfluidity in a Fermi superfluid. The first mechanism is due to 
the excitation of sound modes or phonons. The second mechanism is due to 
the breaking of Cooper pairs, a mechanism which is absent in Bose superfluids. 
These two mechanisms have an opposite behavior as a function of the coupling 
constant, or as a function of the inverse scattering length (as usually is discussed 
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in cold Fermi gases) . While one changes the scattering length across the Feshbach 

resonance ([39||40||4l|[42}|43|[73]) from the BCS limit to the BEC hmit, the binding 

energy of a Cooper pair increases monotonically. On the other hand, the sound 

velocity has an opposite behavior, it decreases monotonically when going form 

the BCS to the BEC limit. As a result, the maximum critical flow velocity 



attains its maximum value approximately at unitarity (76,77), namely Vc ~ 
0.37vf, where vp = fikp/m is the Fermi velocity of a corresponding free Fermi 
gas of the same density. In the following discussion of the vortex dynamics 
Landau's critical velocity will play an important role. In particular, the fact that 
Cooper pairs can and do break up during the dynamics of a fermionic superfluid 
will play also a significant role in our discussion. TDSLDA, unlike traditional 
approaches of superfluid dynamics, naturally allows for Cooper pairs to break up 
when conditions are met. 

We will discuss a class of numerical experiments in which the UFG superfluid 
in its ground state is contained in a vessel with cylindrical symmetry and peri- 
odic boundary conditions along the cylinder axis. The type of containment is 
realized with an almost flat potential and smooth cylindrical walls, high enough 
to prevent the superfluid from ever spilling over. In the first set of experiments 
we introduce adiabatically a "quantum stirrer," parallel to the axis of the cylin- 
der. Such a stirrer can be realized experimentally by a blue detuned laser as 



in the MIT experiment (50). The stirrer rotates with constant angular velocity 
at a radius close to the edge of the system and we vary its velocity from very 
slow to supercritical. After stirring the system for a while the stirrer is extracted 
adiabatically out of the system and the superfluid is left to evolve by itself for a 
while. While the superfluid is stirred energy is pumped into the system mostly as 
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rotational energy stored in newly formed vortices, and part of this excitation en- 
ergy as well as in sound waves. This geometry allows us to simplify the numerics 
by choosing the qpwfs as follows: 

Un{r,t) =^ (i-x.Y>{iknzz)un{x,y,t), Vnif,t) =^ ex.p{iknzz)vn{x,y,t), (6) 

where naturally knz is quantized, due to periodic boundary conditions along the 
cylinder axis. 

If the stirrer is moving very slowly and not enough energy is pumped into the 
system so as to create at least one vortex the superfluid typically returns to its 
initial state at the end of the stirring, see Figure [4j However, above a certain 
velocity, which depends on the geometry of the system one or more vortices 
are created, which tend to form a vortex lattice, constrained by the cylindrical 
geometry of the container however, see Figures [5] and [6| The energy of a single 
vortex state depends logarithmically on the radius of the container, and thus is 
greater in a larger container. Typically sound waves are generated as well, which 
lead to a non-stationary vortex lattice. Surprisingly, a UFG superfluid can remain 
superfluid even if stirred at supercritical velocities, as seen in Figure [6} The 
explanation is relatively simple, even though the occurrence of this phenomena is 
somewhat surprising. A UFG superfluid is a gas, and thus highly compressible. 
During the stirring the gas is gathered by the stirrer as like water by a fast moving 
paddle, and the local density increases. With the density the local critical velocity 
increases as well, and the system can remain superfluid even if it is moving with 
a supercritical velocity of the unperturbed system. When the stirring ends, the 
system comes to an almost steady state, in which the gas is distributed evenly 
throughout the container, mostly along the walls. This is not yet a state in which 
vortices become so close to each other and when a super- vortex state is formed as 
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discussed in Reference (78). The shape of the container plays an important role in 

the nucleation and stability of vortices. From the examples shown on the webpage 



(75), one can see that in an ellipsoidal cylinder the vortex lattice is less stable 
that in a circular cylinder. In particular we could not generate vortex lattice at 
as high stirring velocities in an ellipsoidal cylinder as in a circular one. When the 
superfluid is stirred at velocities well above the critical velocity, as in Figure [7| 
superfluidity is lost, the order parameter monotonically vanishes and the system 
ends up in a normal state. This is an aspect unique to TDSLDA and absent in 
any traditional dynamical models of superfluidity. Naturally, when the system 
becomes normal, the role of dissipative processes is significant and TDSLDA does 
not describe correctly the dynamics of a normal system, and collisions have to 
be included, or a totally different theoretical formalism has to be contemplated 



(65) 



A genuine 3D numerical experiment, illustrated in Figure [8| was performed in 
order to illustrate the generation of vortex rings. An almost impenetrable ball 
was sent flying along the axis of a very long "gun barrel" filled with a superfluid 
UFG at zero temperature. The velocity of the ball was subcritical v = 0.2vc 
as well as subsonic. Example of the dynamics generated by a supersonic ball 
penetrating such a superfluid can be found on the webpage (|75j) . In Figure [s] one 
can see that while traveling the ball is generating perfect vortex rings (mainly due 
to the symmetry of the problem) of various sizes. As in classical hydrodynamics 
( 79p0 ), vortices with larger radii are moving slower than the vortices with smaller 



radii. An attentive reader will also notice that while "colliding" two vortex rings 
do change their radii. 

We have performed a number of numerical experiments to establish microscop- 
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ically in Fermi superfluids the existence of the mechanism suggested by Feynman 



in 1955 (21,81) to be the origin of quantum turbulence. Turbulence in classi- 
cal hydrodynamics exists because there always exist strong dissipative processes. 
Superfluids on the other hand are characterized by negligible viscosity, which in 
the case of UFG even reaches the minimum conjectured value to exist in Nature 
( |71[|70[[72||73p . Feynman conjectured that quantum vortex lines, while moving 
through the superfluid, can often cross and would likely recombine into new vor- 
tex lines. While numerical studies of this mechanism were possible in the case 
of dilute Bose systems, using the Gross-Pitaevskii Equation [T] and classical sim- 



plified models of vortex lines (82,83) or the two- fluid hydrodynamics (22,84), 
that was not the case for Fermi superfluids until now. With the exception of the 
Gross-Pitaevskii equation, in both the vortex filament models and the two-fluid 
hydrodynamics vortex quantization is enforced artificially by hand. The first 
3D microscopic simulations of the quantum vortex nucleation and dynamics in 



a fermionic superfluid were performed in Reference (33), where it was demon- 
strated that 3D quantum vortex lines indeed cross and recombine as Feynman 



envisioned in 1955 (21). We illustrate here the incipient mechanism leading to 
quantum turbulence in a Fermi superfluid with the case of an almost impenetra- 
ble subsonic ball flying parallel to the axis of a very long "gun barrel" filled with 
a UFG at zero temperature. While it traverses the tube, the projectile generates 
a number of imperfect vortex rings and a number of vortex lines with ends on the 
edge of the tube. These vortices propagate at various velocities and often they 
collide with each other and exchange parts of their segments. Sound waves are 
also excited, either by the projectile or during the vortex crossing and recombi- 
nation. In a bigger superfluid volume with a large number of vortex lines these 
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vortex crossings and recombinations lead to quantum turbulence, when the role 

of viscosity is extremely low. 



4.3 Quantum shock waves and domain walls 

Shock waves and soliton waves have been studied for more than a century in 
classical fluid dynamics and plasma physics, and their dynamics and propagation 
properties depend on the presence of dissipative effects and the subtle interplay 
of non-linearites and dispersive effects (85). In the case of UFG moreover, it is 
already an established fact, both theoretically and experimentally, that shear vis- 



cosity is very low, reaching its minimum quantum limiting value (71,70,72,73). 
During the last decade, various experiments performed with dilute cold Bose 
gases and their theoretical interpretation using the Gross-Pitaevskii Equation [T] 
demonstrated that quantum shock waves could be excited in such systems and 



the role of dissipation is negligible ( 86 , 87 , 88 , 89 , 90 ) . In the recent experiment 



performed by Joseph et al. (91) with a cloud of cold ^Li fermionic atoms the exis- 
tence of quantum shock waves in a superfluid Fermi system was confirmed. These 
authors created a long elongated atomic cloud in a harmonic trap. The cloud 
was adiabatically cut in two pieces with a laser beam, which was subsequently 
removed and the two separated clouds accelerated towards each other along the 
long axis of the harmonic atomic trap. Pictures of the collision revealed that af- 
ter merging, two shock wave fronts were formed while the cloud was expanding. 
The experiment was modeled in Reference (^T|) within a modified hydrodynamic 
approach, by adding a phenomenological shear viscosity term 

h + V ■ (vn) = 0, (7) 
nvk + nVk I ^ + /^N + ^1 + i^Vi{n[V v]ki} = 0, (8) 
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where we have suppressed the arguments (r, t). One can choose the mass of the 

atom m = 1, Vk and are the cartezian coordinates of v and V respectively, 

and [V v]ki = Vk^i + Vivi — 25kidiv v/3. Above fi[n] is the chemical potential 

in homogeneous matter at a given number density n, and we have suppressed the 

explicit dependence of the number density n{f,t), velocity v{f,t) and external 

trapping potential U{f, t) on space-time coordinates. Unlike the case of dilute 

Bose gases where the quantum shock waves were interpreted as dispersive shock 

waves ( |86|[87|[88|[89|[90l ) with no need for dissipative effects, the results of the 

experiment ( [9T| ) on UFG received an interpretation similar to classical shock 

waves, when dissipation plays a crucial role in the formation of the shock wave 

front. These two distinct interpretations of the experiments on Bose and Fermi 

dilute gases are difficult to reconcile, especially in the light of what has been 



established so far studying the shear viscosity of these systems (71,70,72,73). 
In the BEC regime, the role of dissipation is negligible and the shock wave and 
the density ripples identified with soliton trains can be described by dispersive 



effects alone. Viscosity was introduced phenomenologically in Ref. (91) to avert 



the onset of a "gradient catastrophe" ( 92 ) . At the same time one would expect 
that in an UFG the role of viscosity is even less important than in a BEC system, 
where viscosity was not needed to model experiments, as the UFG is widely 
accepted as a prime example of an almost perfect fluid. A significant limitation of 
the hydrodynamic approach (91, ) is the inability to describe quantum topological 
excitations (quantized vortices and domain walls in particular), both of which 



have been observed in the similar experiments with bosons ( |86 , |87 , |88||89||90[ ). 

In the case of colliding UFG clouds, we observe the generation of both quantum 
shock waves and domain walls, the excitation of which have been suggested for 
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some time in different kind of simulations ( |93[[94t[95t|96[ ) . The domain walls are 

excitations of the superfluid order parameter and not the number density ripples 

identified as soliton trains trailing the wake of the shock waves, as discussed in 

Refs. ( [86|[87p"8[|89[(90l ) . We will make this distinction in order to avoid confusion. 

We show that the number density of two colliding UFG clouds shows a behavior 

very similar to the one observed in experiment ( j91| . In the wake of the quantum 

shock waves we observe the formation of domain walls. The domain walls emerge 

as quite sharp changes in the phase of the superfluid order parameter by vr, and 

are correlated with minima of the number density. One can distinguish two types 



of domain walls, Fulde-Ferrell-like (97) in which the phase of the order parameter 



changes continuously by vr, and Larkin-Ovchinnikov-like (j98|) when the order 
parameter merely changes signs. Domain walls, which are Fulde-Ferrell-like in 
this case, propagate through the system at slower speeds than the quantum shock 
waves and are topological excitations similar to quantum vortices. Domain walls 
always appear in pairs with opposite jumps of the order parameter phase and 
appear to collide essentially elastically with one another and with the system 
boundary. These phenomena are observed in the absence of any dissipation, 
which is expected to play a negligible role at temperatures close to absolute zero 
and especially in a UFG ([7ll[70}[72}[73]) . 

We have performed simulations of the cold atom cloud collisions assuming that 
qpwfs have the structure as in Equation [6j with periodic boundary conditions in 
the z-direction and a rather stiff harmonic confining potential in the y-direction 
near the box boundary. The time-dependent trapping potential along the collision 



X-axis had a similar profile as the one used in experiment (91 ), namely a shallow 
confining potential in the x-direction with a high potential barrier in the middle 
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that was rather rapidly lowered. The solitons and the shock waves now are two- 
dimensional in character, their stability properties are slightly different than in 



3D (99). Typical results of these simulations are shown in Figures 10 11 12 and 

m 



The simulation results, Fig. 10, show remarkable similarities with the exper- 



iment (91). The shock wave speed in experiment and simulations agree within 
~ 25 — 30%. The differences can be ascribed to various experimental uncer- 
tainties (in particular the particle number and the value of kp after expansion) 
and different geometries. In spite of being confined in the y-direction in a har- 
monic potential, the domain walls are planes perpendicular to the collision x-axis. 



However, in Reference (91) the experimental set-up prevented the authors from 
observing the domain walls. The images corresponding to various frames reported 
there were taken in different realizations of the two colliding clouds. The phase 
differences of the two initially separated condensates are random and cannot be 
controlled from one shot to another, similar to collisions between Bose dilute 



clouds (100). The density profile fluctuations from shot-to-shot in Reference (91) 



point to a rather low spatial resolution attained in these measurements (see Fig. 



2 in Reference (91)), which is insufficient to put in evidence domain walls. We 
have performed simulations by varying the initial relative phase of the conden- 
sates. While the overall picture of the collisions remains unchanged, the number 
of domain walls created varies. The density ripples in the wake of the shock 
waves discussed in experiments with Bose dilute clouds (l86l[87l|88l[89l|90|) and 



interpreted there as a soliton train, are formed here as well. By zooming in the 



online Figures 11 12 and 13 one can notice that before the shock wave is formed 



well defined matter wave interference occurs. The discontinuity in the number 
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density and order parameter at the wake of the shock wave is accompanied by a 



similar discontinuity in the collective flow velocity, see Figures 12 and 13 The 
domain walls which form in the wake of the shock wave have lower speeds. With 
a white circle we marked the region where two domain walls collide, apparently 
elastically. The propagation speed of the domain walls is lower than sound speed. 
One can see the scars left by the domain walls also in the collective flow shown in 



the upper panel of Figure 12 In a hydrodynamic approach, with a phenomeno- 



logical gradient correction term to the EDF, Salasnich and collaborators (101) 



were able to reproduce the shock waves observed in the Duke experiment (91 ) in 
the absence of viscosity. These authors suggest an hydrodynamic approach based 
on the Lagrangean 

?i2|V*(f,t)P 



in-^*{f,t)dt'^{f,t) 



4m 



" ^ 10^ ^ ' ^ -[/(f,t)n(f,t), (9) 

in which they made the identification ^(r, t) = \/n{f, t) /2 exp[ix(r, t)] oc A(f, t) 
and 2mv{f, t) = ?iVx(r, t). Rewritten in terms of n{f, t) and v{f, t) the emerging 
equations look formally like the hydrodynamic equations at zero temperature, 



plus the "quantum pressure" term, which was neglected in Reference (91). In the 
BEC regime this identification amounts to using the Gross-Pitaevskii Equation [T] 
for the dimers/fermion pairs, after replacing the interaction energy as well with 
the dimer-dimer repulsion. In the BCS limit this identification is clearly wrong, 
as well as at unitarity. The hope is that qualitatively one can encapsulate the 
dynamics reasonably well however. These authors also observed the formation 
of density ripples at the front of the shock wave, similar to those observed in 
our simulations, and also in the case of quantum shock waves in Bose systems 
(|86l|87l[88l|89l|90|) . An analytic explanation of the origin of these density ripples 
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in ID Fermi systems was given recently by Bettelheim and Glazman (102). 



4.4 Nuclear collective modes 

The isovector giant dipole resonance (GDR) is perhaps the simplest example of 
a nuclear collective motion of all the protons against all the neutrons. Since its 



observation in the photo-absorption cross section (103), it has been intensively 



studied as it combines several challenging aspects of the physics of the atomic 



nucleus (104,105). Even though the GDR is practically harmonic in character, it 



is not an adiabatic collective mode and various damping mechanisms of the col- 



lective energy are at work (38). In the early models of Migdal (107), Goldhaber 



Teller (108), and Steinwedel- Jensen (jl09') the GDR is described as the relative 



motion of two fluids, either compressible or incompressible, with neutrons and 
protons vibrating around a common center of mass, and the mass dependence of 



the excitation energy reads A~^^^ and A^^^'^ respectively (110). A good estima- 
tion of the GDR vibrational frequency is hu w 80 MeV A~^^^ for spherical nuclei. 
The GDR is interpreted simply as the equivalent of the zero-sound in a nuclear 
system and the size of the nucleus sets a constraint on the largest wavelength. 
In the case of deformed nuclei, the GDR peak is split, with various frequencies 
revealing different principal axes of the nuclear shape. The total width of the 
GDR is mainly due to a couple of mechanisms: the coupling of the GDR to 
complex nuclear configurations r-'-, and the coupling to the continuum, leading 
to the escape of neutrons and protons F^. These two widths contribute to the 
total width of the GDR, F = F-l- + F^, and their relative contributions vary de- 
pending on the mass number A and the N/Z ratio. The escape width is typically 
more important for light nuclei. The physical mechanisms related to F-^ may be 
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quite complicated and involve coupling to low energy surface vibrations, Landau 



damping and collisional damping (106). 

The emerging equations are formally equivalent to the TDHFB approximation 
with local potentials, or to the time dependent Bogoliubov-de Gennes (TDBdG) 
equations: 



ih 



Vkt 
\ J 



A 

-A 



V 



-h* —h* 

''tt u 



h* h* 



Vkt 

\ m J 



(10) 
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-A* 
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where we have suppressed the spatial r and time t coordinates, and k is the label 
of each qpwf [iifco-(r, t), ffcalf) i)]- where a =t,i- The single-particle Hamiltonian 
huf^i{r,t) is a partial differential operator (thus local) and A(r,t) is a pairing 
field, all defined through the normal, anomalous, spin, and isospin densities and 
currents. The interaction with various applied external fields (spin, position 
and/or time-dependent) is described by including the corresponding potentials 
in the single-particle Hamiltonian h„fji[r^t). In one of the first applications of 
TDSLDA to nuclei we have evaluated for the first time the photo-absorption 
cross-section on a trixially deformed open-shell nucleus ^^^Os. Within TDSLDA 
the nuclear response is described by applying an external electromagnetic field 
pulse and performing subsequently a Fourier analysis of the time response in order 



to extract the frequency dependent cross-section ( 37 ) . Typically open shell nuclei 
display superfluidity in both proton and neutron systems. So far there are no 
indications, either experimentally or theoretically, that the two superfluids couple 
to each other in nuclei. The comparison between the theory and experiment 
shows a surprisingly good agreement without any fitting parameters, see Figure 
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14 in spite of the known uncertainties of the nuclear energy density functional. 
It appears that so far the known nuclear energy density functional encapsulate 
reasonably well the gross nuclear properties. 

5 Summary and outlook 

The new framework TDSLDA, which is an extension of the DFT to the real-time 
dynamics of Fermi superfluids incorporates in a natural way the basic elements 
of the hydrodynamic approach at zero temperature (which is nothing else but a 
local implementation of conservation laws), but it also has the necessary quantum 
features to describe a wide range of topological excitations and moreover the 
Cooper pair breaking mechanism. The corresponding EDF is constructed using 
QMC input of homogeneous systems and the properties of the inhomogeneous 
systems are used to validate this functional, or further refine its properties as in 



Reference (52). All required symmetries (translational, rotational, parity, gauge, 
Galilean, isospin, etc.), couplings to various external fields, and renormalizability 
of the theory are naturally incorporated into the theoretical framework. Apart 
from being able to describe correctly known experimental facts, this new approach 
leads to new qualitative predictions (supercritical flow, quantum shock waves and 
domain walls, Higgs modes, vortex crossings, etc.) and it will allow as well to 
study the onset and dynamics of quantum turbulence in Fermi superfluids for 
the first time. It will be particularly interesting to perform further studies of 
the UFG under different conditions, in time-dependent optical lattices with or 
without varying the scattering length both in time and space, and in random 
external potentials, either static or dynamic. It is straightforward as well to 



couple the UFG within TDSLDA to artificial gauge fields (115) and study its 
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response. 

In nuclear physics studies TDSLDA appears as the only theoretical framework 
which would allow us to study on an equal footing both the structure and re- 
actions of medium and heavy nuclei. Traditionally nuclear structure studies are 
performed using different theoretical approaches to describe nuclear structure 
and reactions, and this dichotomy can be put to rest within TDSLDA, in which 
an impinging neutron on a nucleus for example is described on a equal footing 
with the bound neutrons in the target nucleus. In this manner we hope to ad- 
dress for the first time the dynamics of induced fission with neutrons, gamma 
rays, and other projectiles. Properties usually beyond reach within DFT, such as 
various two-body observable, surprisingly, can also be extracted using a method 



suggested quite some time ago by Balian and Veneroni (116117, 118), by imple- 



menting a very simple extension of TDSLDA (65). The dynamics of the vortex 



pinning and de-pinning in neutron star crust (119), and even the onset of tur 



bulent motion due to vortex crossing and recombination (120) are still unsettled 



issues. The vortex pinning mechanism in neutron star crust is still a matter of 



debate (121,122,123) as it is not clear yet whether one has to perform static cal- 
culations of the vortex in the presence of nucleus at constant asymptotic density, 
constant neutron chemical potential or constant particle number in a simulation 
cell. When a vortex is formed a rather large mass redistribution occurs and the 
mass redistribution is further modified in the presence of a nucleus. The uncer- 
tainties arising from using various ensembles in conjunction with uncertainties in 
the nuclear EDF make the resolution of the question if a vortex is pinned or not 
on a nucleus rather unclear. Calculating the actual pinning force requires ana- 
lyzing off-axis configurations that are not stationary. One can characterize the 
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vortex-nucleus dynamics through dynamical simulations, by forcing the nucleus 

to move with an external potential, the vortex (un)pinning can be characterized 



by looking at the rate of energy being transferred to the system (124). 
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Figure 1: The energy of a UFG as a function of the complex pairing field A. 
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Figure 2: The panels o, b and c display response of the homogeneous system to 



various switching time intervals, see Ref. (32) for more details. The initial state 
is the ground state of the UFG, for times t < the coupling constant is slowly 
taken towards the BCS limit in the case shown in panels a and b and at t = the 
unitary value of the coupling constant is restored in a very short time interval. 
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Figure 3: The top two panels display the instantaneous occupation probabili- 
ties of the mode shown in upper panel of Figure [2] corresponding to times t > 
when the pairing field is at its minimum and maximum values respectively by a 
solid blue line with circles. With red dots we plotted the equilibrium occupation 
probabilities corresponding to the same instantaneous values of the pairing gap. 
In lower two panels we show the maximum and minimum values of the oscil- 
lating pairing field and the corresponding excitation energy as a function of the 
frequency of the Higgs-like modes. 
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Figure 4: The magnitude of the pairing field (|A|, top row) and the corresponding 
number density (n, bottom row) for a UFG system composed of 1800 particles 
in a 48^ lattice stirred at OAbVc- No vortices are formed during stirring when the 
system is stirred at a small velocity, however a few sound modes are excited. 
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Figure 5: Several frames showing the UFG at the beginning of the stirring process, 
during the stirring when several vortices were formed, and at the end, after the 
stirrer was extracted and the system was left to evolve by itself. The stirring is 
performed clockwise and the notch close to the boundary of the system shows the 
instantaneous position of the stirrer. Similar situations are illustrated in Figures 
[6]and[7l The magnitude of the pairing field (|A|, top row) and the corresponding 
number density (n, bottom row) for a UFG system composed of 1800 particles in 
a 48^ lattice stirred at a subcritical velocity 0.608vc- Five vortices are formed and 
remain in the system once the stirring concludes. Color coding corresponds to 
blue, lowest, red intermediate and brown the highest value of the corresponding 
quantity. 
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Figure 6: The magnitude of the pairing field (|A|, top row) and the corresponding 
density (n], bottom row) for a UFG system composed of 1800 particles in a 48^ 
lattice stirred at supercritical velocity 1.216fc. Here thirteen vortices are formed 
once the stirring concludes. 
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Figure 7: The magnitude of the pairing field |A| (top row) and the corresponding 
density n (bottom row) for a UFG system composed of 1800 particles in a 48^ 
lattice stirred at supercritical linear velocity 1.824fc- No vortices are formed 
during stirring as the system becomes normal at the end of the run. 
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Figure 8: The magnitude of the pairing field |A| in UFG system composed on 
a 32^ X 96 lattice excited with a centered ball directed in the z-direction. A 
sequence of vortex ring formation is observed. In the sequence of figures time 
progresses left-to-right in a row-major fashion in units of l/ep, and similarly in 
Figure |9j 
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Figure 9: Demonstration of vortex recombination as excited by an off-centered 
ball flying in the z direction. The crossing and recombination of vortex lines were 
conjectured to be the origin of quantum turbulence of superfluids by Feynman in 
1955 dlTl). 
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Figure 10: The left panel displays the two cooling clouds before the collision, 
and the right panel the expanding cloud formed after the collision respectively. 
The expansion front is almost flat, similarly to the one observed in the Duke 



experiment (91) and interpreted as a quantum shock wave. Notice that the 



visual aspect ratio of the clouds is not shown to scale for display purposes. 





Figure 11: Three consecutive frames showing the absolute magnitude of the 
pairing field \ A{x,y,t)\ in the xy-plane at times tep/fi = 30, 350, 690. The x- 
and y-directions (not shown to scale here) have an aspect ratio of ~ 30. We 
have observed the formation of domain walls so far only in traps with elongations 



larger than in experiment (91). 
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Figure 12: The x-component of the collective flow velocity field Vx{x, 0, t) along 
the axis of collision in a space-time diagram. At the front of the two shock 
waves the velocity field undergoes a rapid change in sign, and the matter flows in 
two opposite directions, fully consistent with a shock wave interpretation of the 
expansion front. 
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Figure 13: The left panel shows the phase of the pairing field in units of vr. One 
can clearly see lines of sharp phase change by approximately 27r, consistent with 
the formation of domain walls. The contour of the density profiles is shown in the 
right panel. The white circle shows one area where two domain walls experience 
what appears to be an elastic collision. 
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Figure 14: The calculated photo-absorption cross-section (solid black line), using 
two Skyrme force parametrizations for three deformed open-shell nuclei and the 
experimental (7, ?^) cross-sections (solid purple circles with error bars), extracted 



from Refence (114). With dashed (green), dotted (red) and dotted-dashed (blue) 



lines we display the contribution to the cross-section arising from exciting the 
corresponding nucleus along the long axis, the short axis (multiplied by 2 for 
the prolate nuclei ^''■^Yb and ^^^U) and the third middle axis in the case of the 
triaxial nucleus ^^^Os. 



